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Abstract 

The impact of an isolated vortex in a compressible Keplerian disk is examined using higher order 
numerical solutions of the Euler and entropy-conserving Energy equations. The vortex is 
stretched by the background shear flow with longer lasting anticyclonic vortices persisting for 
about 10 vortex revolutions. Simultaneously, the vortex emits transient radial waves consisting 
mainly of ax i symmetrical weak shock waves and a slower, nonaxisymmetric Rossby wave. 

These waves may contribute to certain transient events in protoplanetary disks. The vortex 
stretching and waves were found to have little long-term feedback on the baseline “standard solar 
nebula” disk structure and confirm the extremely stable structure of non self-gravitating disks. 


I Introduction 

Light, cool disks with approximately Keplerian rotation curves are common models for 
protoplanetary disks. Such gaseous disks form the background for solid body accretion and gas 
dynamic processes that no doubt plays a critical role in the planet formation process 
(Weidenschilling & Cuzzi 1993, Beckwith et al. 2000). Earlier work on idealized infinitely thin 
disks (Lynden-Bell & Pringle 1974, Adams et al. 1989) has evolved to investigations of multi 
dimensional disks with vertical and radial mean flow structure. (Klahr et al. 1999) There is now 
an extensive literature base concerning the dynamical properties of such disks. Furthermore, a 
number of disk-related wave systems were identified as being relevant to the solar nebula. A 
survey of the field to the mid 1990s by Papaloizou & Lin 1995 emphasizes the role of the 
underlying gas and seriously consider wave motion as one option for mass and momentum 
transfer. 

The strong synergy between geophysical and astrophysical problems was recognized decades 
ago. In particular, coherent vortices as an important flow structure in protostellar disks was 
advocated by v. Weizacker 1945 and examined in some detail by Adams & Watkins 1995. 
Bracco et al. 1998, 1999 consider transport properties of vortices in barotropic fluids using 
spectral-based finite-difference techniques with direct application to protoplanetary disks. They 
showed that coherent vortices last long enough to form “lumpy structures” that may concentrate 
solid objects and accelerate planetesimal formation. Recently, Godon & Livio 1999 extended the 
work to a compressible viscous nebula by simulating discrete vortices in a viscous accretion disk 
using related pseudo-spectral numerical methods. Using a polytropic gas model they find that 
anticyclonic vortices last from 10 to 100 revolutions depending on the chosen viscosity. No 
acoustic or shock wave motion is reported, but their equations clearly support such waves. Barge 
& Sommeria 1995 investigate long-lived vortices in a disk and predict that particles may migrate 
to the center of a vortex after a few vortex revolutions so the issue of vortex dynamics in the 
particle accumulation process is of current interest. This paper will consider vorticies in gaseous 
disks and the dual role of the underlying disk medium in propagating a variety of wave systems 
throughout the disk and in simultaneously stretching the vortex by localized velocity shear. 
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Sheehan et al 1999 examine the analogy between meteorological and cosmological wave 
motions in some detail, emphasizing the role of Rossby waves. They argue that vortices and 
Rossby waves probably play a role in the protoplanetary nebula. Davis et al 2000 (hereafter 
DSC) further investigate the role of vortices and Rossby waves in a model solar nebula They 
consider the evolution of a compact region of coherent vorticity using an incompressible flow 
model and the vorticity transport equations. They found that compact vortices induce weak 
Rossby waves that slowly radiate over the entire disk. The vortices ultimately dissipate, but 
anticyconic vortices were shown to be more persistent than cyclonic vortices. 

In this paper we expand the analysis of vortices to compressible disks and show how a vortex, 
once formed, can generate disturbances that propagate to the outer reaches of the disk. We use a 
standard model of the solar nebula with a variable thickness and power-law radial variations of 
the equilibrium variables. We show that discrete vortices evolve in two directions: first into a 
ring structure at the radial location of the vortex and second into two separate wave packets. One 
is a fast wave based on the acoustic time scale and the second a slower Rossby wave motion. The 
fast wave evolves into a weak cylindrical shock that envelops the entire disk. We find that no 
instability mechanisms are activated and the waves have little impact on the underlying disk 
structure once they traverse the disk. 

A third order in space/ fourth order in time explicit numerical algorithm is developed to solve the 
“column averaged” Euler and Energy equations. The equations are cast in a semi-conservative 
form designed to exploit the underlying eigenvalue structure. The base flow assumes mechanical 
equilibrium with assumed density, pressure, and local disk height based on standard solar nebula 
models. A seeded compact vortex is chosen whose circulation is approximately the rotation rate 
of Jupiter and diameter appropriate for the thin disk approximation. 

With respect to the question of vortex generation in a protostellar disk Li et al 2000 derive 
dispersion relations for unstable Rossby wave phenomena induced by a strong local entropy 
(temperature) gradient. Li et al. 2001 simulates this instability and subsequent anticyclonic 
vortex formation using numerical solutions of the Euler and energy equations. The standard solar 
nebula model used in this paper does not include mean variations in entropy and, as a result, only 
globally stable wave systems were found. 

In section II the numerical problem is formulated based on a thin disk flow model that is initially 
in mechanical equilibrium. A new derivation of the three-dimensional equilibrium flow and the 
definition of the column-averaged quantities are presented in an Appendix. A semi-conservative 
form of the Euler and Energy equations is described and used to develop a high-resolution flux- 
biased finite difference algorithm. In Section m the vortex evolution is examined using 
instantaneous and azimuthally averaged flow quantities. Finally, the relevance of this model to 
some open questions regarding disk evolution is discussed. 


II Problem Formulation 

(a) The equilibrium flow 
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The equilibrium surface density and pressure from given distributions on the central-plane are 
derived in the Appendix. These column-averaged thermodynamic quantities define a fictitious 
polytrope (Hunter, 1972) that is related to the physical gas by the parameters in Table Al. The 
temperature of the disk is conveniently defined from the column-averaged perfect-gas equation 
of state in terms of the column pressure P and the surface density a: P/Po = (ctT)/(GoTo). The 
initial central-plane thermodynamic quantities po. To, n and y are all given constants. Here we 
will take values representative of a “standard solar nebula.” Let ro be 1 AU and n = -3 as 
recommended by Weidenschilling 1977. Cuzzi, et al, 1993 used po= 1.4xl0‘ 9 g/cm 3 and To = 280 
K at ro= l. From the perfect gas equation of state (assuming a molecular weight of 2.34 g/mole), 
po = 13.83 dynes/cm 2 . The specific heat ratio y is taken as 1.5 and the value of GM appropriate 
for the solar nebula. 

It is useful compute the dependent variables in astronomical units as the compressible flow 
equations contain a number of thermodynamic quantities that are not conveniently 
nondimensionalized. In addition, the key role played by the sound speed is expressed in a more 
physical context. Length, mass, and time are taken as AU, solar mass (SM), and years (yr). 
Certain parameters become extremely simple in this system. For example, G = 4 tt and M = 1 . 
Table 1 indicates the major equilibrium flow parameters at several radial locations using the data 
in the Appendix. 

TABLE 1 


Equilibrium Flow Parameters 


r, AU 

H, AU 

a, SM/AU 2 

P, SM/yr 2 

T, K 

a, AU/yr 

1 

0.082 

2.05E-04 

7.73E-06 

180 

0.274 

5 

0.409 

4.11E-05 

3.09E-07 

36 

0.123 

10 

0.817 

2.05E-05 

7.73E-08 

18 

0.087 

15 

1.226 

1.37E-05 

3.44E-08 

12 

0.071 

20 

1.634 

1.03E-05 

1.93E-08 

9 

0.061 


The second column is the half height of the nebula using formulas in the Appendix. The 
computed wave motions should have wavelengths larger that these values to avoid three- 
dimensional effects. The surface density, pressure, and temperature decrease quite rapidly. The 
last column is the local sound speed which decays in this particular case as r , the same decay 
rate as the basic Keplerian flow so the azimuthal Mach number stays constant at about 23. The 
entropy is conveniently measured by the polytropic “constant” K’ = P o‘ r It is initially a global 
constant of about 183 in astronomical units. The ratio of gravitational (GM/r) to thermal (P/o) 
energy, which governs the disk thickness, is also constant and approximately equal to 1050. 


(b) Equations of motion and solution method 
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The equations of motion can be written in a variety of ways. Analytically based asymptotic 
methods were very popular in the 1950s and 60s and used conventional non-conservative 
equation sets. In the continuing search for effective numerical methods, the conservative form of 
the equations seems preferable. Unfortunately, astrophysical problems include gravitational 
body-force terms that upset the symmetry of the classical conservation equations. In addition, 
when expressed in non-Cartesian coordinates the fully conservative equations become quite 
complex. Davis 2001 compares straightforward use of the classical conservation equations in 
curvilinear coordinates and showed that the effect of rotation could be significant. The numerical 
error when computing a simple Keplerian flow in polar coordinates was examined and traced to 
an incompatible finite-differencing sequence. 

The alternate form of the Euler/Energy equation is shown below. It is based on the fact that 
astrophysical flow fields are driven by differences between two large quantities (gravitational 
and centrifugal forces). If these relatively large quantities can be expressed as source terms, 
incompatible discretization errors will not occur. If u r and uq are polar velocities, e t the total 
energy, p the pressure, and p the density, then 
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where the dependent variable Q = [p, pu e , pu e , e,] and the second and third vectors are fluxes. 
This equation has two important properties. First, a baseline Keplerian flow is satisfied exactly 
by the undifferentiated source terms. Second, the (r, 0) eigenvalues can be shown to be [u r , u r , 
u r +c, u r -c] and [ue, ue, ue+c, ue-c] where c is the local sound speed. This makes it easier to 
construct wave-following algorithms in direct analogy with Cartesian systems. The physical 
nature of the solar nebula is such that perturbation radial velocities are close to sonic while the 
baseline azimuthal velocity is at quite a high Mach number (ue/c ~20). This form of the Euler 
equations has been presented previously (Mair et al 1988; Godon 1997) but the favorable 
properties mentioned above seem not to have been utilized. 

The equations are integrated in time using an explicit numerical algorithm. Higher order methods 
are necessary when computing evolving waves on reasonably sized grids. A critical feature is the 
need to transmit waves with minimal phase error (dispersion). The change in Q over a small time 
increment is computed using a fourth-order Runge-Kutta method. Spatial derivatives are 
computed with a third-order biased upwind scheme (a five-point stencil with three upwind and 
one downwind mesh point). The spatial stencil was applied to a flux-split form of eq (1). The 
two flux vectors can each be decomposed into the product of a 4x4 “Jacobian matrix” and the 
vector of dependent variables Q. The Jacobian matrix, in turn, can be converted by a similarity 
transformation into a product of three 4x4 matrices consisting of a pre-matrix, a diagonal matrix, 
and a post-matrix. The elements of the diagonal matrix corresponding to 0 and r are 
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[u e ,u g ,Ug + c,Ug-c ] and [u r ,u r ,u r +c,u r -c] respectively. These are the allowable propagation 

speeds of the waves at each point in space-time. In the former case u e will dominate c since the 
azimuthal flow is highly supersonic. The radial velocities are approximately sonic so they will be 
more sensitive to the ambient sound speed. The sign of each of these speeds is examined and 
those of like sign are combined in each coordinate direction. The matrices are then recombined 
into two flux vectors that reduces to the matrix equation: 


a _ a 

r — 0-\ 

dt dd + 


F + + 


d0- 


■ F H G* h — ~G + R — 0 

3r + dr 


where derivatives with superscript + (-) are backward (forward) biased. The most 
computationally extensive part of the algorithm is computing the flux vectors which is required 
at each stage of the Runge-Kutta time advancement algorithm (four times per time step for each 
coordinate direction). This finite difference scheme is contrasted with other conservation-law- 
based schemes aimed at precise shock wave resolution. They are usually not of high-order 
accuracy away from shock waves which is the trade-off for highly resolved shock waves. The 
development of higher order algorithms with precise shock resolution is the subject of intense 
current research. 


Appropriate boundary conditions for this class of problems is formidable. These highly sheared 
flows propagate waves whose wavelength is comparable with the shear rate and commonly used 
radiation boundary conditions are not applicable. Considering the lack of reliable and robust 
numerical radiation boundary conditions, a second order boundary extrapolation was used. This 
simple approach was appropriate for this problem. Wave reflection from the inner boundary was 
not a problem. (However, disturbances seem to focus on the boundary, which is probably a 
numerical artifact.). The outer boundary was taken far enough away so that reflections did not 
yet interfere with the outwardly moving wave. Periodic boundary conditions were imposed in the 
azimuthal direction. 


Mechanical equilibrium was disturbed by a vortex placed at 4 AU. Vortex parameters were 
chosen to satisfy the thin disk approximation. The vortex is distorted by the local equilibrium 
flow and simultaneously transforms some of its energy into an outwardly propagating trailing 
spiral. The computational domain consists of a disk ranging from 1 to 24 AU using (0, r) = 
(65,300) mesh points. Additional computations were made on (125,300) and (65, 600) point 
grids with no significant changes to the flow. Both cyclonic and anticyclonic vortices were used 
as initial flow perturbations. 


Ill Analysis of Vortex-induced Waves 

The imposed vortex is a transient event and induces an excess velocity field upon the slightly 
modified mean Keplerian flow of eq (A5). This excess flow is shown in Figure 1 along with the 
initial patch of vorticity (maximum value about -1 .0 yr 1 ) in the third quadrant. Its core diameter 
is about 1 AU and peak azimuthal velocity of about .14 AU/yr. This corresponds to a circulation 
(T = 27tr V onVevort) about 50 times that of Jupiter which represents a relatively strong event in the 
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protoplanetary nebula. The local sound speed, and disk thickness at this radius are about 3.14 
AU/yr and 0.3 AU respectively. The vorticity is quite compact, but note the large extent of the 
vortex-induced flow including significant radial velocities. The vortex imposes a thermodynamic 
(pressure, density, and temperature) response on the disk. The entropy (s = P <3' r ) is initially 
constant everywhere except for a small region near the vortex. Subsequently, the flow evolves by 
redistributing the initial circular vortex into a sheared vorticity patch and a variety of propagating 
waves. Unlike the incompressible disk considered by DSC, the compressible disk supports a 
much more varied and complex flow structure. 

Early (a few vortex revolutions) and later (many vortex revolutions) flow fields are expected to 
be associated with different physical effects. The former event is a relatively rapid relaxation and 
redistribution of the vortex caused by the dominating Keplerian shear at 4 AU. Wave radiation 
processes are expected to dominate later. A very important parameter relating to wave systems is 
the smaller flow-induced radial velocity. In the remainder of the paper we will consider both 
effects. These flow events are illustrated using contour diagrams and instantaneous flow patterns 
in sec HI(a) and with azimuthally averaged quantities in sec 111(b)- 

(a) Instantaneous Flow Patterns 

Figure 2 compares perturbed vorticity contours (scaled on the instantaneous peak vorticity) for 
cyclonic (initial vorticity aligned with the Keplerian vorticity) and anticyclonic vorticity. The 
horizontal scale represents angular position (0 < 0 < 360) and the vertical scale is the physical 
radius (3 < r, AU < 10). The first pane] in each part of Figure 2 corresponds to the vortex in 
Figure 1. The cyclonic vortex in Figure 2(a) shears away quite rapidly (note contours at the 
radius of the initial vortex) and simultaneously emits a vorticity wave that propagates radially at 
the local sound speed. This axisymmetric wave is still in the field of view at 8 revs but by 12 
revs it has moved beyond 10 AU. The initial vortex has almost completely dissipated by 12 revs. 
In contrast, the anticyclonic vortex maintains it coherence for a much longer time span (Figure 
12(b)). The vortex is still coherent, but it also emits a vorticity wave into the outer disk. The 
contours in Figure 2 are drawn relative to the peak vorticity at each time (much larger in the 
latter case) so the propagating waves are actually of comparable amplitudes to those in part (a). 
This is clarified in Figure 3 with instantaneous radial vorticity profiles at 4 revs. The peak 
anticyclonic vorticity is larger than the cyclonic vorticity. The propagating waves, on the other 
hand, have similar amplitudes. These results show that the compressible disk responds in two 
ways; first a local flow subject to intense velocity and secondly a wave field that is relatively 
insensitive to source details. 

Compressible and incompressible simulations are compared in Figure 4. Part (a) shows contour 
maps for an anticyclonic incompressible vortex from DSC with similar characteristics to that of 
Figure 2(b). The same contour lines are used and the stretched vorticity near the source relative 
to the compressible flow is apparent. (The incompressible Rossby waves studied in DSC are very 
low level and do not appear with this choice of contours.) Even after 4 revolutions the 
incompressible vortex has stretched well beyond that of its compressible counterpart indicating 
that vortices in compressible flows retain coherence for longer times. The difference between the 
extended vorticity fields in incompressible and compressible flow can be explained by the fact 
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that vorticity is not a conserved quantity in the latter case. However, the potential vorticity 
(called generalized vortensity by Adams & Watkins 1995), defined by the ratio of vorticity to 
surface density is a conserved quantity. Figure 4(b) shows contours of the compressible potential 
vorticity using the same progression of contour values. It is much more compact and is a closer 
analog to the incompressible vorticity. The role of potential vorticity in particle accumulation in 
the protoplanetary nebula is a subject for future study. 

The evolution of the wave radiation processes is shown in the following figures. Figure 5 shows 
vorticity and surface density contours superimposed on perturbed velocity vector after 1.33 
revolutions. The perturbation vorticity has stretched significantly from its initial value and is now 
in the fourth quadrant in part (a) of the figure. During this early phase the flow is rapidly and 
continuously sheared by the baseline Keplerian flow. A positive induced radial velocity (note the 
velocity vectors on the outermost spiral) indicates an outward moving wave. The trailing spiral 
moves clockwise and outward under the influence of the baseline Keplerian flow and the induced 
velocities. The leading edge of the spiral is in the first quadrant at this instant. A coherent 
overdensity (which grew from an initially insignificant value) is visible in part (b) at the same 
location as the vorticity. At the same time, there is a strong density wave forming in the second 
and third quadrant. This strong density wave is the result of nonlinear distortion of the transient 
pressure and is the first indication of shock wave formation. 

Readjustment of the density wave from an isolated disturbance to an axisymmetric shock wave is 
shown in the density contours of Figure 6. The sequence illustrates the tendency of a supersonic 
Keplerian shear flow to '‘circularize” isolated disturbances. The flow evolves as a system of 
trailing spirals that move very rapidly in the azimuthal direction (close to Keplerian speed, above 
Mach 20) and with a radial speed close to the local sound speed (Mach 1). The outermost spiral 
is a precursor compression wave (light spirals in the outer disk) followed by a pair of shock 
waves. The precursor evolves radially near the local sound speed, but the shock wave moves 
faster since it is continuously overtaking slower moving material. Beyond the shock wave the 
low level precursor wave is starting to be reflected by the imposed boundary conditions at 24 
AU. Behind the shock wave is a slower wave system with a nonsymmetrical amplitude 
distribution. This wave system is identified as a Rossby wave based on its similarity (generation 
time and radial wave speed) with those computed in DSC. Of these three waves, the main 
disturbance is confined to the pair of shock waves steadily moving outward. By 12 revolutions 
the wave system has enveloped almost the whole disk. A perspective view of the entire density 
field at this time is shown in Figure 7. The two shock waves are readily apparent along with the 
slower Rossby waves. There is also a residual density collar surrounding at the initial vortex 
radius. Space-time features of these waves are further analyzed into component azimuthal 
modes. 


(b) Averaged Flow Variables 

A finite Fourier series representation is used to define azimuthal averages. For example, the 
mean surface density (m = 0 mode) and fundamental (m = 1 mode) are: 
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where N is the number of azimuthal grid points and 9j = 2Ji(j-l)/N. Nonlinear evolution of the m 
= 0 mode is the main pathway for radial transmission of information while the m=l mode with a 
cos 0 dependence is useful for characterizing the slower Rossby wave component. The 
axisymmetric perturbations of surface density and entropy characterize the wave field. The 
quantities Co(r, t) and s 0 (r, t) are shown in Figures 8 and 9 at selected times to 100 yrs. These 
waves characterize the fastest signals emitted by the deforming vortex. The time scale for 
information transfer is directly related to the sound speed; it takes about 12 revolutions of the 
vortex at 4 AU to affect the region out to 20 AU. 

The surface density in Figure 8 evolves from essentially no disturbance to a propagating density 
wave in about 10 yr. The evolving wave consists of a precursor compression/expansion followed 
by two sharp compressions and a decaying wave packet. The precursor wave decay is probably 
governed by geometrical spreading in the cylindrical disk. The compressions evolve into shock 
waves that maintain their amplitude. The stronger shock separates a weak leading shock and the 
trailing wave packet. This discrete wave system is the only long-range signal emitted by the 
continuously decaying vortex. A persistent nonlinear density remnant remains at 4 AU. 
Differences between rapid compressions and shock waves are not easily discemable from density 
profiles, but entropy changes are sensitive indicators. Figure 9 clearly shows the formation, 
growth, and decay of the shock waves as measured by the entropy. The first jump appears after 
20 yr. at 6 AU, grows to 60 yr., and then begins to decay. Meanwhile the second stronger shock 
forms behind the first and feeds on its predecessor. The entropy signal is slightly contaminated 
by a numerical artifact from the finite difference equations but does not affect the physical nature 
of the evolving flow. The leading shock move at a speed ranging from .14 AU/yr at 38.4 yrs to 
.10 AU/yr at 96 yrs, both values being greater than the local mean sound speed. The two shock 
waves at 96 yrs possess temperature jumps of .003 and .873 deg K respectively, a fraction of the 
ambient temperature of approximately 12 deg K. Also note steadily falling entropy in the region 
of 4-6 AU. This is where the slower Rossby waves appear and an azimuthal average smears out 
the discrete nature of these disturbances. 

The slower moving Rossby waves are most visible in the potential vorticity, the compressible 
counterpart of the conserved vorticity. Figure 10 shows a space-time graph of the fundamental 
(m =1) component of this quantity. The shock wave positions are clearly marked along with a 
variable wavelength disturbance that lags the faster moving acoustically generated waves. The 
slow waves decay rapidly and are similar in nature to the Rossby waves in DSC. 

The computed wave systems are summarized in Figure 1 1 against a background of characteristic 
curves belonging to the undisturbed flow. The family of dotted characteristics correspond to the 
trajectories dr/dt = a(r) = const with their slopes proportional to the inverse sound speed. The two 
solid lines estimate the shock wave path obtained by fitting a parabola to the density peaks at 
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each indicated time. The shock speed decreases due to a combination of its propagating into an 
inhomogeneous medium and nonlinear wave interactions. Their slope is always less than the 
sound speed indicating a mild local supersonic shock speed. At larger radii the mean vorticity 
rapidly decays and acoustically generated waves should be indistinguishable from classical 
cylindrical weak shock waves. Considered as an inverse signaling problem, the nature of the 
source (as a vortex) cannot be reconstructed from this far field generic signal. The dashed lines 
in the figure trace maxima of the fundamental mode potential vorticity (a representative trace 
from Figure 10 is redrawn for clarity showing the location of both waves). Their slope is always 
greater than that of the background characteristics indicating a much slower phase speed. These 
waves remain in the vicinity of the source disturbance about twice as long as the shocks, but 
represent a weaker wave system. 


Conclusions 

Adams and Watkins 1995 present a comprehensive and still relevant overview of the many open 
issues related to vortices in protoplanetary disks. Here we address a few of the issues using a 
numerical model of the nonlinear evolution of vortices in disks. The main findings relate to 
bifurcation of the vortex structure into propagating wave packets and into an axisymmetric ring 
structure. In particular the potential vorticity (e.g. vorticity per unit surface mass) emerges as a 
useful structure in its own right. 

The main conclusions from this investigation are: (1) As a result of Keplerian shearing motion a 
vortex persists for about 10 revolutions. Whether or not this is sufficient time for a critical mass 
of particles to accumulate needs further investigation. Even though vortex persistence is 
dependent on a number of parameters, it is unlikely that vortices of this size could last more than 
decades of orbits. However, it is clear that anticyclonic vortices are definitely more robust. 

(2) Another outcome of this investigation is the non-participation of the vortex in the baseline 
structure of the disk. Consistent with classical notions regarding wave propagation in non 
participating media, the wave moves through the disk without materially affecting it. Global 
mass accretion and angular momentum transfer in protostellar disks are associated with large- 
scale, widely-spaced asymmetrical spirals of matter while the wave packets computed here are 
restricted to tightly-wound spirals and axisymmetric shocks. 

(3) The response of disks with other equilibrium profiles may be quite different. If self- 
gravitating effects are included (heavy disks, e.g. Laughlin et al 1998) vortex-induced 
instabilities could “light up” the disk with strong transient signals and drive accretion and 
angular momentum transfer processes via gravitational torques. Another possibility is that 
thermal instabilities due to thickness gradients or radial entropy gradients could induce 
significant transient events. In any event it is hard to see how any of these time-limited effects 
could sustain global turbulence in the disk and be active for long epochs. The only possibility 
would be an almost continuous barrage of external vortex-inducing events. 
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Appendix 

Mechanical Equilibrium 
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A polytropic gas rotates only under the influence of a central gravitator. An approximate 
equilibrium solution is presented in Ruden et al. 1988 and Boss 1993. Here we will describe the 
mechanical equilibrium as a general solution of the governing partial differential equation. In 
spherical coordinates the three non-vanishing flow quantities are the rotational velocity V<p, the 
density p and the pressure p. All three are independent of the azimuthal coordinated cp. The 
system is closed by a polytropic equation of state p(R,0) =K p(R,0) Y where K is the polytropic 
constant and y is the ratio of specific heats. If the pressure is eliminated, the two remaining 
equations that conserve rotational and radial momentum are: 


■ pV^ cot 6 + Kyp 7 1 = 0 


dd 


V, 


'L +K7p r-*£ + 5M£-= 0 


R 


dR R 


(Al) 


Eliminate the rotational velocity to obtain a first order partial differential equation for the 
quantity p 7 " 1 : 


r - 1 


R 1 cos 6 


dp_ 
dR 


Y — 1 dp 7 

— — GM cos d - R sin d — ^ — = 0 
Ky dd 


(A2) 


The general solution of this partial differential equation for p 7 1 involves an arbitrary function of 
the quantity R sin0 that is just r, the radial distance from the axis of rotation. 


P 


Y - 1 


f(R sin0)[l- 


7-1 CM ( 1 J_ 

Ky f(Rsin0 ) Rsin# R 


(A3) 


Identify f(r) with the (y- 1 )st-power of the arbitrary central-plane density p c (r) to complete the 
solution. The pressure is easily obtained from the polytropic law and the azimuthal velocity from 
the first of eq (Al). The final relation among the three quantities are: 


GM d\ogr 


|i/(M> 


yKr- 0 


vhR,e) = 


GM 


Rsind 


(A4) 
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where ac, is the isothermal sound speed defined by the square root of p c /p c - A reasonable 
approximation to eq (A4) in cylindrical coordinates is to let R ~ r and expand 1/r-l/R in a Taylor 
series about z = 0. The final mechanical equilibrium state for a thin disk is: 


y-\GM/rz 

p = pm) i -h- — i — 

2 Y “a r 


-il/(y-l) 


= P<(r) 




1- 


y-lGM/rz 2 
2 Y *1 


Y?(Y- 1> 


{ | 70 a r ^jOg Pc 

GM d log r 


(A5) 


The coefficient of the term multiplying z 2 /r 2 proportional to the ratio of the local gravitational 
potential energy GM/r to the thermal energy squared. For a typical protoplanetary nebula this 
ratio is O(10 3 ) so the thermodynamic quantities decrease very rapidly from z = 0. The azimuthal 
velocity, on the other hand, depends only on the radius r and is independent of z, which is a 
general statement for a barotropic fluid (Lebovitz 1979, p 236). The density and pressure decay 
rapidly from the midplane while the rotational velocity remains very close to the Keplerian value 
throughout the flow field. (The coefficient of the corresponding term in is now the inverse of 
the energy ratio mentioned above.) These formulas are the same as those used by other 
investigators concerned with vertical equilibrium in polytropic disks. 

The “column-averaged” density (a) and pressure (P) are obtained by integrating the quantities 
over a disk whose half-height H is the inverse of the coefficient of z 2 in the above formulas. The 
equations for the surface density and surface pressure are (Hunter, 1972, p 224): 




-H 

H 


H 2 


]/(r-n 


p = p, Mj 


-H L 


1 — 


H 2 


rur- 1 ) 


dz = Hp c {r)^[\-x 2 ] Kr 0 dx = fiHp c {r) 


(A6) 


-I 


where a and (3 are defined by Gamma functions that depend only on the ratio of specific heats y. 


Now let the variation of the central density and pressure be simple power laws of the radial 
coordinate; p c (r) = p 0 (r/r 0 ) m and p c (r) = p 0 (r/r 0 ) n . It follows that the variation of H obeys the 
power law H(r) = Ho (r/ro) M where Ho = (2y po/po /(y-l)/GM) 1/2 and |1 = (3-m+n)/2. The surface 
density and pressure now have spatial variations of the form <T(r) = apo Ho (r/r 0 ) m+M and P(r) = 
Ppo Ho(r/ro) n+ ^. This is a parametric representation relating surface density and pressure in terms 
of r. Eliminating r implies the relation between a and P: 
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p= & . (j m+ f = K'cr r 

n+p 

{(XH 0 )m+^ 

which indicates that the surface thermodynamic quantities are also polytropes with a modified K 
and a different specific heat ratio. The exponents m and n are constrained so that m y - n — 0. 

In summary, mechanical equilibrium is defined by the given quantities po , po? to, n, GM, and y. 
The derived quantities are m, ot, (3, p, T, and Ho- Table AI shows representative values of those 
parameters in terms of y for the case when n = -3. 

TABLE A 1 


Values of the Equilibrium Parameters (n = -3) 


Y 

m 

a 

P 

11 

r 

1.2 

-2.500 

0.739 

0.682 

1.250 

1.4 

1.3 

-2.308 

0.876 

0.786 

1.154 

1.6 

1.4 

-2.143 

0.982 

0.859 

1.071 

1.8 

1.5 

-2.000 

1.067 

0.914 

1.000 

2 

1.6 

-1.875 

1.137 

0.958 

0.938 

2.2 
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List of Figures 

1 . Vortex location and induced velocity field at initial instant of time. The vortex is located at 
about 135 deg on the r = 4 AU circle. 

2. Contours of equal strength vorticity relative to minimum (maximum) vorticity at each time 
step. Time is measured by vortex revolutions at 4 AU. (a) Initial cyclonic vortex, 
perturbation vorticity > 0. (a) Initial anticyclonic vortex, perturbation vorticity < 0. 

3. Vorticity profile along a ray at 163 deg showing decay of initial vorticity at 4 AU and 
instantaneous location of the wave field. Shock waves have already formed in the region 6-8 
AU. 

4. Contours of (a) vorticity in an incompressible disk model, (b) potential vorticity in a 
compressible disk model. Anticyclonic vorticity shown using the same scale as Figure 2. 

5. Contours of (a) vorticity and (b) surface density after 1.37 revolutions of the initial vortex. 
Rapid shearing and initial wave systems are apparent. 

6. Evolution of the surface density over the entire simulated disk, 1 < r, AU < 24. (a) Initial 
flow, (b) 4 revs, (c) 8 revs, (d) 12 revs. 

7. Perspective view of the total surface density at 12 revs. (Note: perturbed density drawn to 
twice the scale of the equilibrium density for clarity.) 

8. Profiles of mean perturbation surface density Go(r, t) at selected times showing global shock 
wave propagation in the disk. 

9. Values of mean perturbation entropy so(r, t) indicating discontinuities at the shock wave 
position. 

10. Potential vorticity profiles for the m = 1 mode showing cos 0 dependence of slower moving 
Rossby waves. 

11. Space-time diagram of shock wave trajectories (solid), Rossby wave tracks (dashed), and 
equilibrium flow characteristic curves (dotted). 
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